02. 比对 |
您所在的位置:网站首页 › samtools flagstat › 02. 比对 |
比对
这一步将使用hisat2或STAR对数据进行比对。hisat2在比对速度与资源占用上说比STAR要好上不少,但是也有不少人建议使用STAR,因为STAR可以获得更好的比对结果。因此具体的选择还是根据实际情况定。关于hisat2、STAR、topHat的对比可以看这篇文章。 HISAT2 分步进行原则上,hisat2先将数据比对到sam文件,再使用samtools转换为bam文件后再进行排序。其中-x后跟的是索引文件的名字(不包括扩展名)。hisat2对SRR1374921比对用时3分钟,比对率在98.24%。 1 2 3 4hisat2 -p 10 \ -x genome/mm39 \ -U rRNA/SRR1374921_filtered.fq.gz -S bam/SRR1374921.sam再使用samtools转为bam 1 2samtools view -bS bam/SRR1374921.sam \ | samtools sort -@ 10 - -o bam/SRR1374921.bam最后再用samtools进行排序并创建索引 1samtools index bam/SRR1374921.bam然后其实可把SRR1374921.sam删掉了。 合并进行其实我们也可以一条命令,直接生成排序后的bam。使用管道和tee同时生成比对报告。 1 2 3 4 5 6 7hisat2 -p 10 \ -x genome/mm39 \ -U rRNA/SRR1374921_filtered.fq.gz \ | tee >(samtools flagstat - > bam/SRR1374921.hisat2.flagstat.txt) \ | samtools sort -O BAM \ | tee bam/SRR1374921.bam \ | samtools index - bam/SRR1374921.bam.bai写一个循环来跑这些数据。下载的内容里有个mergelist.txt,有所有数据样本名,把样本名截下来。 1cut mergelist.txt -d'_' -f1 > samplelist.txt然后循环运行 1 2 3 4 5 6 7cat samplelist.txt | while read line; do \ hisat2 -p 10 -x genome/mm39 \ -U rRNA/${line}_filtered.fq.gz \ | tee >(samtools flagstat - > bam/${line}.hisat2.flagstat.txt) \ | samtools sort -O BAM \ | tee bam/${line}.bam \ | samtools index - bam/${line}.bam.bai; done比对完成后,进入下一步统计reads数。 STAR同样的数据,STAR需要50min才能完成,比对率85.76%。 1 2 3 4 5 6 7 8STAR --genomeDir genome/star_index/ \ --readFilesIn rRNA/SRR1374921_filtered.fq.gz \ --runThreadN 10 --outSAMtype BAM SortedByCoordinate \ --quantMode GeneCounts \ --outFileNamePrefix bam/SRR1374921. \ --readFilesCommand zcat mv bam/SRR1374921.Aligned.sortedByCoord.out.bam bam/SRR1374921.star.bam samtools index bam/SRR1374921.star.bam对文件存放结构不满意的需要把log移动到log文件夹。 由于在索引中还加入了GTF,STAR会同时输出比对到各个基因的counts,结果在ReadsPerGene.out.tab中,后三列分别是“counts for unstranded RNA-seq”, “counts for the 1st read strand aligned with RNA (htseq-count option -s yes)”, “counts for the 2nd read strand aligned with RNA (htseq-count option -s reverse)”。根据需求选用哪一列。 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |